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ABSTRACT 

Cyclic spectral analysis is a signal processing technique designed to deal with 
stochastic signals whose statistics vary periodically with time. Pulsar radio emission 
is a textbook example of this signal class, known as cy do stationary signals. In this 
paper, we discuss the application of cyclic spectral analysis methods to pulsar data, 
and compare the results with the traditional hlterbank approaches used for almost 
all pulsar observations to date. In contrast to standard methods, the cyclic spectrum 
preserves phase information of the radio signal. This feature allows us to determine 
the impulse response of the interstellar medium and the instrinsic, unscattered pulse 
profile directly from a single observation. We illustrate these new analysis techniques 
using real data from an observation of the millisecond pulsar B1937+21. 

Key words: methods: data analysis — pulsars: general — pulsars: individual (PSR 
B1937+21) — ISM: general — scattering 



1 INTRODUCTION 

Pulsar radio emission produces some of the most 
information-rich signals found in astronomy. Variation ex- 
ists over a wide range of time scales, from ns-duration gi- 
ant pulses (|Hankins et al.ll200ot ) , to tiny decade-scale timing 
fluctuations that may eventually r esult in the first di rect de- 
tection of gravitational radiation (|Hobbs et al.ll201oh . Prop- 
agation of the radio pulses through the ionized interstellar 
medium (ISM) adds another layer of wavelength-dependent 
complexity (and information), including dispersion, multi- 
path scatterin g, and diffractive and refractive scintillation 
|Rickettlll99(il ). These effects provide a unique tool for ex- 
ploring the structure of the ionized ISM in our galaxy on 
a wide range of scales. However, they also represent a still 
not fully understood source of s ystematic error for high- 
precision pulsar timin g projects jHemberger fc Stinebring 
l20Qg| ; IColes et al.ll20ld ). 

The complexity of the pulsar signal has led to the inven- 
tion of a variety of specialized signal processing techniques. 
Currently, the most common approach to pulsar signal anal- 
ysis is to divide the digitally-sampled baseband voltage into 
a number of frequency channels using a fast Fourier trans- 
form (FFT) or autocorrelation. The signal power is then 
detected in each channel, and integrated over a short time 
to produce a power spectrum estimate. This results in well- 
known tradeoffs and limitations on time and frequency res- 
olution due to the competing eff ects of dispersive smearin g 
and inverse channel bandwidth (|Lorimer fc Kramer|[20(jl . 
A major breakthrou gh came with the introd uction of co- 
herent dedispersion l|Hankins fc Ricketdll975T ). in which a 



filter is applied within each channel pre-detection to com- 
pletely remove the effect of dispersion, for a known disper- 
sion measure (DM). This allows wider channels and hence 
higher time resolution to be achieved, an d is commonly used 
for high-precisi on pulsar timing (e.g., I Hot an et alj 120061 : 
lDemorestll2007j ). This method is still bound by the funda- 
mental time/frequency resolution relationship AtAv > 1, 
with the values of At and Av fixed at the time of observa- 
tion. 

Traditional spectral estimation procedures are based on 
an implicit assumption that the input signal is stationary, 
or at least approximately so over the spectrum integration 
timescale. However, the pulsar radio signal can be more ac- 
curately described as cyclostationary - it is a non-stationary 
random si gnal whose sta tistical properties vary periodically 
with time (Rickett 1975). Over the past 20 years, statistical 
analysis techniques have been developed specifically to deal 
with cyclostationary signals, the study of which is known 
as cyclic spe ctral analysis (see reviews by iGardnen Il99ll : 
lAntoni|[2007l ). These methods examine not only power ver- 
sus frequency as in conventional spectral analysis, but also 
the correlation between spectral components that arises in 
cyclostationary signals. Cyclostationarity has recently been 
explored in a radio astronomy context for th e identificatio n 
and removal of radio frequency interference (Fcliachi 2010). 
However, these methods have not previously been applied 
to the pulsar signal itself. Doing so provides access to sig- 
nal content that is unreachable via standard methods, and 
opens the door to a variety of novel analysis techniques. 

In this paper, we describe the theory and practice of 
applying cyclic spectral analysis techniques to pulsar data. 
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Section [2] gives a brief overview of cyclic spectral analysis 
and introduces the primary quantity of interest, the cyclic 
spectrum and its Fourier transforms. In Section [3] we de- 
scribe in detail how the cyclic spectrum can be computed 
for pulsar observations, and present an example using ac- 
tual data. Finally, Section [4] shows how the pulsar cyclic 
spectrum can be used to measure and correct for ISM scat- 
tering, and discusses possible future applications of these 
methods. 



2 CYCLOSTATIONARITY AND CYCLIC 
SPECTRAL ANALYSIS 

We begin with a brief overview of stochastic processes and 
spectral analysis, to establish terminology. This is followed 
by a introduction of the main features of cyclic spe ctral anl- 
ysis. T he more rigorous treatmen ts presen t ed by IPapoulisl 
for stochastic processes and lGardnerl J 1993 ) for cyclic 
spectra are useful references. While we frame the discussion 
in terms of continuous-time functions, the discretization of 
the following equations is straightforward. 

2.1 Standard Spectral Analysis 

A general stochastic signal, x(t), can largely be characterized 
by its second-order statistics. In the time domain, these are 
correlations, defined here symmetrically: 

C x (t,T) = E{x(t+^)x*(t-^)} (1) 

Here, the expectation value represents an average over many 
realizations of the signal. Stationary processes are those 
whose correlation function depends only on the time dif- 
ference r between any two points and not explicitly on the 
time t. While the signals themselves have random variation 
versus time, their statistics are constant. The signal can al- 
ternately be characterized by the power spectrum S x (v), 
which is the Fourier transform of the correlation function 
C x {t). 

For ergodic signals, the correlation function can be es- 
timated from a single realization of the signal via a time 
average: 

a(r) = T~ 1 j\(t+ T -)x*{t- T -)dt (2) 

This forms the basis for practical spectral analysis of ran- 
dom signals. In practice, both ergodicity and stationarity 
are almost always assumed a priori, since one usually has 
only one realization of a signal. Due to the Fourier-pair un- 
certainty relations, the spectral resolution Av is limited to 
approximately T _1 or larger. 

In radio astronomy, x(t) is the baseband voltage, a 
quantity proportional to the incident radio wave's electric 
field, mixed (frequency-shifted) to near zero frequency. Early 
digital radio spectrometers often explicitly evaluated Eqn. [2] 
using dedicated hardware. In modern systems, S x (v) is usu- 
ally computed using a fast Fourier transform (FFT) of x(t). 
In order for a pulse to be resolved in time, T must be small 
compared to the pulse width. For pulsars, this puts a limit on 
the finest achievable frequency resolution of Av > (PA<jf>) _1 
where P is the pulse period and A(f> is the pulse width in 
turns. 



2.2 Cyclic Spectral Analysis 

An interesting class of nonstationary signal are those whose 
statistics vary periodically with time. These are known as 
cyclostationary signals and their correlations exhibit period- 
icity: 

C x (t, r) = C x (t +P,t) = C x (t mod P, r) (3) 

It is important to note that although C' x (t, r) is periodic 
in t, it is not generally periodic in r, and x(t) itself is not 
necessarily a periodic function. Simple examples of cyclosta- 
tionary signals can be obtained by periodically amplitude- or 
frequency-modulating a stationary noise process. In this pa- 
per we will sometimes use the quantity pulse phase <f> = t/P 
in place of t as an independent variable in expressions such 
as C x (c/>,t). 

The cyclic analog of the power spectrum can be ob- 
tained by Fourier transforming C x (t,r) along both the t 
and r axes: 

+ oo P 

S x {u;a) = P" 1 j drfdt C x (t, T ) e -™ VT e-' M ° ct (4) 

-oo 

This quantity is known as the cyclic spectrum and it is a 
function of two frequency variables. The frequency v is con- 
jugate to r exactly as in the standard power spectrum. The 
cycle frequency a is conjugate to t, and takes on the discrete 
values ct n = n/P - it is computed here as a Fourier series 
rather than a continuous Fourier transform. In the station- 
ary signal case, S x (v; a) = for a ^ 0. 

The periodic correlation can be estimated from a sin- 
gle signal realization similarly to Eqn. [2] by "folding" the 
correlations modulo the periodicity P: 

N 

C x (t,r) = N- 1 x(t + nP + I)x'(t + nP - I) (5) 

This again assumes ergodicity of the signal, over N peri- 
ods. In this computation, the integration time T is given by 
NP, the number of periods folded. The integration time sets 
an ultimate limit on the range of r that can be measured, 
and again the spectral resolution is limited to Av > T^ 1 . 
However, in constrast with Eqn. [2] T is not constrained by 
the pulse width, as the nonstationarity has already been cor- 
rectly taken into account in Eqn. [5] This decouples the pulse 
phase resolution from the frequency resolution. 

It is easy to see that various symmetry relationships 
exist among these quantities. In our symmetric formulation, 
C x (t,-r) = C x (t,r). Similarly, S x (v;-a) = S%{v;a). A 
useful intermediate quantity is the periodic spectrum S x (v, t) 
that is obtained by Fourier transforming C x (t,r) with re- 
spect to r alone. By symmetry, S x (v, t) is purely real- valued, 
and therefore is straightforward to visualize. S x (v,t) is su- 
perficially similar to signal power as a function of time and 
frequency. However as discussed below, it also contains in- 
formation about the signal phase, and thus is not strictly 
positive. 

It can be shown that the cyclic spectr um also repre- 
sents correlati ons in the frequency domain (|Gardnerlll99ll ; 
lGardnerlll992h . Using the Fourier transform of x(t), X(y), 
we have: 

S m (u;a) = E{X{y + a/2)X*{v-a/2)} (6) 
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This relationship directly leads to a key result on which 
we will rely heavily in later sections: If y(t) is the result 
of passing x(t) through a linear, time-invariant filter with 
impulse response h(t) (frequency response H{v)), 

y(t) = h(t)*x(t) (7) 

Y{u) = H{u)X{u) (8) 
then the cyclic spectra of x and y are related by: 

S v {v;a)=H{v + a/2)H*{v-a/2)S x {v;a) (9) 

Eqn. [9] shows that the phase content of H (v) is pre- 
served in the cyclic spectrum, whereas in a standard power 
spectrum the only information retained is the filter magni- 
tude \H{y)\\ This is a fundamental difference between the 
two approaches that will enable most of the novel applica- 
tions described in £[4] 

In practice, computation of any of the cyclic spectrum 
variants will be done using band-limited, digitized sampled 
values. That is, the signal from the radio frequency of inter- 
est will be filtered with a bandpass filter of width B, mixed 
to baseband and critically sampled at the Nyquist rate. It 
is well known that any remaining power outside the base- 
band frequency range —B/2 < v < B/2 will be aliased 
into the sampled values. Similar considerations applied to 
cyclic spectra show that the resulting S(v; a) is only valid 
within a diamond-shaped region around the origin, defined 
by \a/2\ + M < B/2. 



3 IMPLEMENTATION FOR PULSAR DATA 

A simple but useful model of the pulsar radio signal y(t) as it 
is rec eived on Earth is amplitude-modulated noise |Rickettl 

y(t\ to, v ) = hisAi(t; to) * (p(tfi(£o); v )rn(t)) (10) 

Here, nt (t) is a stationary noise process representing the in- 
trinsic pulsar radio emission. p(<f>; vo) is a periodic function 
describing the average pulse profile shape at a radio fre- 
quency vo- The timing model Q(to) gives the apparent rota- 
tional frequency of the pulsar at a given time. This includes 
the slowly-changing intrinsic spin frequency of the pulsar as 
well as Doppler and relativistic terms due to the motions 
of the Earth and pulsar. The impulse response function, 
hisM(t;to), describes the effect of the ionized interstellar 
medium (ISM) on the radio signal. It contains both disper- 
sion and scintillation/scattering terms, and slowly evolves 
with a characteristic timescale td due to the relative mo- 
tions of the Earth, ISM, and pulsar. For clarity, we will of- 
ten omit the slow dependences on to and vo in the following 
discussion. 

Under the approximation that the average intrinsic pul- 
sar emission spectrum S n (v) = So is flat (valid over small 
fractional bandwidths), the pre-ISM pulsar cyclic spectrum 
takes the simple form S x {y; a n ) = I(n)So, where I(n) is the 
Fourier transform of the intensity profile I{4>) = [p(<^)] 2 and 

1 For simplicity, the model presented here neglects the complex 
"subpulse" behavior of many radio pulsars as well as other fil- 
tering processes caused by the ionosphere, troposphere, telescope 
equipment, etc. This does not alter our basic conclusions. 



a n = nQ. The addition of the ISM filtering process produces 
an expression for the final pulsar cyclic spectrum: 

S y {u;a n ) = H ISM {v+ ^)H* ISM (v - ^)I(n)S (11) 

The duration of the ISM impulse response, Tism, can be 
approximated as the combination of the two primary effects, 
dispersion and scattering: Tj SM ~ t'dm + T sc- ln order to 
resolve the full frequency structure of Hism{v), the cyclic 
spectra must be computed with a frequency resolution Ais < 
(2tism)~ ■ Since Hism(v) is slowly evolving with time, the 
cyclic spectrum integration time should be kept shorter than 
the diffractive scintillation timescale td for the signal model 
presented in Eqn. II II to apply. 

3.1 Time-domain Approach 

A time domain approach to computing the cyclic spectrum 
is based on directly evaluating Eqn. [5] In this method, the 
complex-conjugated baseband values are delayed by an in- 
teger number of samples ("lags"), and then multiplied by 
the original undelayed version. The pulse phase for each 
point in this cross-multiplied data stream is calculated us- 
ing standard methods (polynomial expansion of Q(to)) and 
then they are folded into the desired number of pulse phase 
bins. The process is repeated for the desired number of de- 
lays (Ni a g). Due to symmetry, only positive lags need be 
evaluated. When combined with the "mirrored" negative lag 
values, the result is an Nbin by N c han = 2Ni ag — 1 array rep- 
resentation of C x (4>,t). This computation can be preceded 
by standard coherent dedispersion at a known DM to reduce 
the required frequency resolution, effectively setting tdm to 
~ 0. 

This time domain method has the benefit of straightfor- 
ward handling of the relationship between the slowly chang- 
ing pulse period and sampling rate. It also is easily inte- 
grated into standard pulsar observing practices: The "zero- 
lag" term is simply a standard folded pulse profile as would 
usually be computed by a traditional pulsar processing sys- 
tem. The drawback to this method is that the computational 
cost per sample scales linearly with N c h an , quickly becoming 
impractical if fine frequency resolution is desired. For real 
time computation over a total bandwidth B with frequency 
resolution Af, the computational cost in floating-point op- 
erations per second (flops) is given by 

Ct(B,Au) ~4N pol — (12) 

Here, N po i equals 2 if only total power terms are computed, 
or 4 if polarization cross-products are included. An imple- 
mentation of the time-domain cyclic spectrum method is 
now included in the open-so urce DSPSR software packag^3 
(|van Straten fc Bailesf[201oh . 

3.2 Frequency-domain Approach 

It is also possible to compute the cyclic spectrum in the fre- 
quency domain, based on Eqn. [S] In this method, the raw 
samples are first Fourier transformed with resolution deter- 
mined by the tism timescale. The frequency domain data 

2 http://dspsr.sourceforge.net 
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Figure 1. Periodic spectrum S x {v,t) of B1937+21 using fre- 
quency resolution of 4 MHz/6320 = 0.642 kHz and pulse phase 
resolution of 1.55 ms/511 = 3 fis. The boxed region is shown in 
detail in Figure [2] 



is multiplied by an appropriate phase gradient to align the 
data in pulse phase, and a coherent dedispersion filter is 
optionally applied at this point. The frequency domain cor- 
relations of Eqn. [6] are then computed for N[ ag "spectral 
lags" of spacing a n . The number of spectral lags determines 
the pulse phase resolution, with Nbin = 2N' lag — 1. However, 
since there is generally not an integer number of pulse peri- 
ods per FFT length, the frequencies required for evaluating 
Eqn. [6] fall between FFT bins. In order to avoid artifacts in 
the resulting cyclic spectra, the FFT output must be inter- 
polated to obtain the correct fractional frequencies. Here we 
use a 7-point sine interpolation with good results and leave 
a detailed study of the effect of this interpolation for future 
work. Finally the correlations are averaged into N c han fre- 
quency bins. As in the time domain case, applying coherent 
dedispersion reduces the final required frequency resolution. 

The frequency domain method has the advantage of 
the computational requirements scaling linearly with Nun 
rather than N c han as in the time domain version. This is at- 
tractive for situations where very fine frequency resolution 
is desired. However, there is some added complexity due 
to the frequency interpolation, and this method is a much 
larger departure from standard pulsar observing strategies. 
For this method, the number of flops required for real-time 
computation is 

C F {B, Av) ~ N bln (2N mt + AN pol )B + 5Blog 2 (13) 

where Nmt is the number of points in the interpolation ker- 
nel. The second term is the approximate cost of the initial 
FFT. Real-world FFT performance can depart significantly 
from this scaling depending on the implementation details 
and computing hardware architecture. 



3.3 Hybrid Filterbank Methods 

Neither the time-domain method of Section 13.11 or the 
frequency-domain method of Section 13.21 is efficient for di- 
rect application to a wide (B > 10 MHz) frequency band. Ct 
scales as B 2 , and while Cf formally grows only as BlogB, 
the FFT length required for dedispersion grows as B 2 and 
eventually becomes impractical in current computing hard- 
ware. The computational burden can be reduced by first 
splitting the wide band into a number of sub-bands via a fil- 
terbank, and then applying either the time or frequency do- 
main cyclic spectrum computation method in each sub-band 
separately. This is similar to how most coherent dedisper- 
sion algorithms currently work. However, in contrast with 
standard filterbanks, to compute a continuous cyclic spec- 
trum over the full wide band, the sub-bands must be made 
to overlap in frequency. As explained at the end of Section 
12.21 the cyclic spectrum can not be properly computed for 
frequencies within a/2 of an edge of a band-limited signal. 
Therefore to avoid gaps at each sub-band edge, the sub- 
bands must overlap by an amount B a i a max , and the in- 
valid edges discarded before finally appending the sub-band 
results together. This approach results in a total number of 
real-time flops given by 

C H ~Cfb(B,B s ) + — ?-—Ct.f(Bs, Av) (14) 

where Cfb is the cost of the filterbank operation, B 3 is 
the full width of each sub-band (including the overlapped 
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Figure 2. Left: Detail of periodic spectrum of B1937+21. The black rectangles each contain unity time-frequency area. The fixed 
resolution implied by any one of these, as is produced by a traditional filterbank, would miss various features in the data. Note also the 
appearance of both positive (dark) and negative (light) signal regions at fine scales. This is due to the phase conten t of Hjsm(u) present 
in th is representation. Right: The same data processed with a standard coherently dedispersed filterbank (DSPSR; van Strat enfc Bailesl 
l2010t) , using a frequency resolution of 12.5 kHz, and corresponding intrinsic time resolution of 80 fis (black rectangle). The data are 
oversampled in pulse phase by a factor of ~20, resulting in the same number of profile bins as in the left panel. In both panels the 
zero-point has been determined by measuring the mean off-pulse level, and the data have been normalized by the peak value. 



portion), and Ct,f is the cost of either of the previously 
explained methods applied to a single sub-band. Due to 
the relationship B s > B a i ^ a max , the minimum possible 
sub-band width is set by the desired pulse phase resolution 
A0= {2Pa max )- 1 . 

3.4 Example Data 

We illustrate the concepts just discussed with an example 
of cyclic spectrum computation on actual pulsar data. We 
observed the 1.55-ms pulsar B1937+21 at Arecibo Obser- 
vatorjjf] usin g the Astronom y Signal Processor (ASP) pul- 
sar backend (|Demorestll2007l ) . Dual-polarization, 8-bit base- 
band voltage samples were recorded in a 4 MHz band cen- 
tered on 428 MHz. Cyclic spectra were computed via the 
frequency domain methocfl discussed in £13.21 using a FFT 
length of 542288 points, 256 spectral lags and integrated for 
Ti n t ~ 2 minutes. The FFT length was chosen based on the 
dispersion time tdm ~ 30 ms. A coherent dedispersion filter 
using a DM of 71.019 pc cm -3 was applied before summing 
the cyclic spectra to the final output frequency resolution 
of N cha „ = 6230 (Au = 0.642 kHz). When transformed to 
the periodic spectrum domain (Figures [T] and [2}, the 256 

3 The Arecibo Observatory is part of the National Astronomy 
and Ionosphere Center, which is operated by Cornell University 
under a cooperative agreement with the National Science Foun- 
dation. 

4 The same data analyzed with the time domain method pro- 
duced identical results. 



harmonics result in Nbi n = 511 [PAcf> = 3.0 ps). The spec- 
tra were computed separately for each polarization, then 
summed to make the total intensity spectrum shown. For 
this example, the radio frequency resolution Av was chosen 
to be approximately equal to 1/P, providing equal spacing 
between points in the v and a dimensions in the cyclic spec- 
trum. 

PSR B1937+21 is an ideal case for showing the utility 
of these methods. Near 430 MHz, the typical ISM scatter- 
ing ti me is r se ~ 40 fis, comparable to the intrinsic pulse 
width (|Cordes et al.lfl990T ). Even if coherent dedispersion is 
first performed, the typical filterbank approach forces one to 
choose between: 1. Resolving the scintillation structure (on 
frequency scales ~ t^, 1 ) but losing almost all pulse shape 
information; or 2. Retaining enough time resolution to re- 
solve the pulse, but missing most of the fine frequency struc- 
ture. These tradeoffs are graphically illustrated in Figure [2l 
where it can be seen that the cyclic spectrum approach re- 
tains information on both small pulse phase and frequency 
scales. As previously mentioned, the periodic spectrum con- 
tains signal phase information and the on-pulse region can 
be seen to contain both positive and negative components, 
relative to the mean off-pulse level. Also notable are the fine 
"streaks" that appear at later pulse phases. One possible ex- 
planation for these features is that in this part of the profile 
the signal is dominated by more highly delayed components, 
making the periodic correlation C(<f>, r) contain signal power 
at larger |r| values than are present at earlier phases. When 
transformed to the periodic spectrum domain, this results 
in finer frequency structure appearing at later pulse phases. 
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4 APPLICATIONS 

Given the signal model presented in Eqn. llll it is possible to 
determine both the ISM response and intrinsic pulse profile 
directly from a single cyclic spectrum. The two-dimensional 
cyclic spectrum S{v; a n ) contains N c h an x N[ ag data values, 
while Hism{v) and I(n) are described by only N c h an + N[ ag 
model parameters. This provides sufficient constraints for 
both to be determined via iterative least-squares minimiza- 
tion. A detailed analysis of this met hod will be presented in a 
separate paper l|Walker et al.ll201ll ). One prev ious measure- 
ment method for hisM(t) has been published (|Walker et al.l 
2008), based on a dynamic spectrum phase retrieval pro- 
cedure. The cyclic spectrum method is much simpler since 
in this case the wave phases can be measured directly, giv- 
ing an estimate of hiSM from a single "snapshot" obser- 
vation. The cyclic spectrum also incorporates pulse profile 
shape information, which is lost in standard dynamic spec- 
tra. This naturally leads, for the first time, to a true coher- 
ently descattered pulse profile shape (Figure |3). In contrast 
with previous inte nsity-based profile deconvolution methods 
(|Bhat et al.ll2004T l. this requires no assumption of a specific 
functional form for Hism and is not affected by ambigu- 
ity between intrinsic and ISM-induced profile features. The 
ISM response shown in Figure has an initial exponential 
decay followed by a more slowly-decaying tail. This will be 
interesting to compare in detail with the pre dictions of the 
stand ard Kolmogorov scattering model (e.g., iRickett et al.l 
l2009h . 

There are several degeneracies that the descattering 
process alone can not resolve. As is clear from Eqn. 1111 mul- 
tiplying Hism by a constant phase factor will produce no 
change in the observed cyclic spectrum. Similarly, without 
additional assumptions, the pulsar's intrinsic flux So is de- 
generate with the magnitude of Hism- Most critically for 
pulsar timing, an arbitrary rotation can be applied to I((j>), 
and absorbed into Hism- That is, at this level of analysis 
it is impossible to distinguish an ISM-induced delay from a 
pulsar spin deviation or other timing effect. Resolving this 
situation to obtain properly scattering-corrected timing will 
require the development of additional analysis techniques. 
This could range from assuming a constrained form or ap- 
plying moment analysis to hisM to physical models of the 
spatial distribution of scattering material, and is an active 
topic for further study. Compared with previous methods, 
the cyclic spectrum provides a qualitatively new way to mea- 
sure the ISM response, and this new information dramati- 
cally expands the possibilities for scattering corrections to 
timing data. 

In addition to the determination of Hism, Eqn. fallows 
for the applcation of other phase-coherent filters to the final 
integrated cyclic spectra. This technique can be used to per- 
form coherent dispersion corrections, within the limits of the 
cyclic spectrum resolution. The high frequency resolution 
of the pulsar cyclic spectrum enables precise post-detection 
removal of narrow-band radio frequency interference, with- 
out sacrificing pulse phase resolution. More sophisticated 
approaches may also incorporate infor mation gained from 
the cyclostationarity of the interference (|Feliachi||2010T ). 

Although the discussion in §S\ focused on the analysis 
of a single stochastic signal, pairs of correlated cyclostation- 
ary signals can be analyzed via cyclic cross-spectra, in corn- 
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Figure 3. Top: Absolute square of the impulse response 
\hisM{~t)\ 2 determined from the data shown in Figure [T] Bot- 
tom: Corresponding standard (dashed) and scattering-corrected 
(solid) pulse profiles. All functions have been normalized by their 
peak value. The x-axis range covers the same total time span in 
both panels. 

plete analogy with standard cross-spectra (|GardnerHl992l ). 
For dual-polarization radio data, this results in the creation 
of cyclic Stokes parameters, and allo ws for the applica tion 
of phase-coherent matrix convolution l|van Straterj|2002l ) di- 
rectly to the cyclic spectra. For radio interferometers, ana- 
lyzing data from antenna pairs will produce cyclic visibil- 
ities. Along with recently developed VLBI imaging tech- 
nique s for investigating pulsar scintillation (Briske n et al.l 
2010), this could prove to be an extremely powerful tool for 
exploring the ISM. 



5 CONCLUSIONS 

Cyclic spectral analysis is a powerful new observational tech- 
nique for studying radio pulsars. It provides a data represen- 
tation that simultaneously preserves both high pulse phase 
resolution and high frequency resolution information about 
the signal. With the accompanying preservation of signal 
phase content, this allows fundamentally new analysis tech- 
niques not possible with standard filterbank data. This holds 
promise both for increasing our understanding of the ion- 
ized ISM, and eventually for removing ISM scattering as 
an obstacle to achieving the highest possible pulsar timing 
precision. 
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